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Abstract 



H. 

C/3 , Modified estimators for the contribution rates of population eigen- 

^ values are given under an elliptically contoured distribution. These 

■ estimators decrease the bias of the classical estimator, i.e. the sample 
contribution rates. The improvement of the modified estimators over 
the classical estimator is proved theoretically in view of their risks. We 

£f~^ ■ also checked numerically that the drawback of the classical estimator, 

■ namely the underestimation of the dimension in principal component 
analysis or factor analysis, are corrected in the modification. 
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1 Introduction 

Let U be the population covariance matrix of a p-variate random vector 
). Let A = (Ai, . . . , A p ), Ai > • • • > X p denote the eigenvalues 
of U, then the population contribution rates are defined as 

v 

r = (n, . . . ,t p ), T,i = Xi/^Xj, i = l,..., p. (1) 

i=i 

The contribution rates play an important role in statistical linear models. 
Especially in principal component analysis or factor analysis, it gives an 
important information for the determination of the model's dimension; "How 
many principal components substantially represents the total variance ?" is 
a basic quercitin in principal component analysis, and the number of factors 
to be incorporated in a model is a crucial problem in factor analysis. For 
this issue, the most simple and widely used methods are the following ones 
based on the population contribution rates. 
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1. The cumulative percentage of the eigenvalues 

With a cut-off t*, we determine the smallest integer m for which 



m 



(2) 



i=i 



to be the number of principal components or factors to be retained. 
Practically a number between 0.7 and 0.9 is often chosen as a cut-off 



2. The relative size of each eigenvalue 

If the ith eigenvalue is larger than the average of the population vari- 
ance X]f=i \/P> ^ ne corresponding principal component or factor is to 
be retained. This criteria is equivalent to check whether r, satisfies the 
inequality 



This is also equivalent to "Kaiser criterion" in factor analysis, which 
asserts that the number of the eigenvalues larger than unit of the pop- 
ulation correlation matrix should be the number of factors. 

Many methods have been proposed for the determination of dimension re- 
lating to principal component analysis or factor analysis (or more generally 
covariance structure model). See Jolliffe (2002) and Fabrigar et al.(1999), 
both of which give an extensive review of the methods for choosing a dimen- 
sion respectively for principal component analysis and exploratory factor 
analysis. According to Jolliffe (2002)'s classification, there are several cate- 
gories for the methods other than that based on the population contribution 
rates; 1) Hypothesis testing method, 2) Information theoretic method , 3) 
computer-intensive method. 

Furthermore it might be better to add another category, "large dimen- 
sional random matrix method", if we could name it. For the past decade, 
while the results have (re) accumulated on "general asymptotics" that con- 
siders the limiting operation of both p (the dimension) and n (the sample 
numbers), we have seen much improvement on this method. The limiting 
distribution of the sample eigenvalues under general asymptotics gives some 
novel ideas for the dimension determination. It is appealing that several 
simulations show that the arguments based on general asymptotics is effec- 
tive even if n and p are relatively small. See e.g. Kritchman and Nadler 
(2004), Ulfarsson and Solo (2008). They propose highly efficient methods 
for a so-called "spiked covariance model" (see the equation (13"8l ). For large 
dimensional random matrix theories used in these papers, see the references 
therein. We also refer to Paul (2007), Nadler (2008), Karoui (2009) for more 
recent developments. 



t*. 



n > P 



-i 



(3) 
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We should notice that the concept of "dimension" could be rather ambigu- 
ous term. In the fields such as physics or chemistry, it is often the case that 
there exist "signals (components)" and "noise" in its own mechanism with 
clear distinction. Naturally the analysis of the covariance structure is aimed 
at "detecting" the numbers of the signals (components), as the term "signal 
detection" indicates. On the contrary in psychology or economics, a "factor" 
in its theory is rather abstract object and sometimes impossible to draw the 
line between the "factors" and "noise" . There we could only say some factors 
are trivial while the others are nontrivial. Hence the number of the factors 
(i.e. dimension) are not considered to preexist before the statistical infer- 
ence but rather be determined through the inference so that we can carry 
out dimension reduction without serious loss of information. We could say 
it is "deciding" the dimension. Considering the both cases, "detecting" and 
"deciding" the dimension, it seems that there is no single excellent method 
that is universally effective. After all we had better choose effective methods 
according to the purpose of the determination of the dimension and/or the 
presupposed mechanism of data generation. 

Back to the methods (J2J), (EJ) of our concern, their cut-off values seem 
somewhat ad-hoc without rigorous theoretical background. We naturally 
raise a question such as "Why is 0.9 for t* chosen ?" We only could insist 
that it is nothing more than conventional criteria for the inference like a 
given significant level (e.g. 5%) in a hypothesis test. Nevertheless, they 
have been widely used and incorporated into many softwares for statistical 
analysis because of their simplicity and easiness for calculation. They do not 
suppose any rigid data generation mechanism, which is often suitable for the 
purpose of dimension "decision (reduction)". We think that the improved 
inference on the population contribution rates could make some contribution 
to the better dimension reduction. In this paper we focus ourselves to the 
point estimation of r using the sample covariance matrix. 

Let A denote the (unbiased) sample covariance matrix and l\ > ■ ■ ■ > 
I* > be its ordered eigenvalues. Then the sample contribution rates are 
defined as 



Traditionally (and perhaps almost always) the set of sample contribution 
rates has been used for the estimation of the population contribution rates. 
Hereafter the sample contribution rates as an estimator of r will be called 
"classical estimator" and denoted by = (f[°\ . . . ,fp ^), where ff = 
di, i = l,...,p. 

As far as we know, for the estimation of r there is no other option than the 
classical estimator. However seems to have nonnegligible bias. It is well 
known that the sample eigenvalues I* = (1*, ■ ■ ■ ,1*) are biased. Specifically 



p 




1 < % < p. 



(4) 
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saying, Aj, i — 1, . . . ,p is majorized by E(l*), i — 1, . . . ,p, that is, 

m m 

< 5>(J?), l<m<p. 

3=1 3=1 

This fact makes us conjecture that is also biased. Let S = (N — 1)A, 
where N is the number of the samples. In the case S is distributed as 
Wishart matrix and A's have no multiplicity, the expected value of di is 
expanded with respect to the degree of freedom n{— N — 1) as follows (see 
the proof in Appendix); 

E((]) T , -J ^g^jj 2 A? , Aj v A, 1 
\ (Y7j=i x j) 3 (E^A^ + E^A^At-A,/ 

+ 0(n- 2 ) 

(5) 

The coefficient of the n~ l term is complicated, but we easily notice that when 
Aj's are close to each other, large positive (negative) bias might take place for 
the smaller (larger) i's. Note that the similar expansion with respect to both 
p and n might be possible in view of "general asymptotics" . Please refer to 
Nadler (2008) for a matrix perturbation approach. 

The simulated results of the distribution of d^s under the condition U 
is a identity matrix can be found in Mandel (1972) and Krzanowski (1979). 
They observe the large bias of d^s. (See also Sugiyama and Tong (1976), 
Konishi (1977), and Huang and Tseng (1992) for the distribution of d^s.) 
Table 1 shows the simulated values of E(d i ) ) i = 1, . . . , 10 calculated from 
10000 random 10-dimensional Wishart matrices with the degree of freedom 
30 generated under several patterns of A. (The total sum of Aj's always 
equals one, hence A« = Tj, 1 < % < p. ) As the figures in the table show, 
it is not seldom that that the bias of the first or last few £ , (d i )'s surpasses 
50% (sometimes 100%) of Aj, while the sign of the bias for the middle part 
of E(d, i ) , s is quite unstable. 

The aim of this paper is to derive an alternative estimator that modifies 
the bias of the classical estimator. In the next section, first we show the 
distribution of the sample contribution rates is identical under a class of 
elliptical distributions. Second we propose a class of new estimators and 
show their superiority to the classical estimator under the class of elliptical 
distributions from a decision theoretic point of view. In the third section, by 
simulation studies, we clarify other preferable aspects of the new estimator. 

2 Main Result 

2.1 Framework 

Before deriving new estimators, we formulate the estimation problem of our 
concern. Let x^\i = 1, . . . ,N be independently and identically distributed 
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Table 1: Bias of dj 



Ai 


A 2 


A 3 


A 4 


A 5 


A 6 


A 7 


A 8 


A 9 


Aio 


E(di) 


E(d 2 ) 


E(d 3 ) 


E(d A ) 


E(d 5 ) 


7~ 1 / T \ 

E(d 6 ) 


7~ 1 / T \ 

E(d 7 ) 


7~ 1 / T \ 

E(d 8 ) 


7~ 1 / T \ 

E(d 9 ) 


7~ 1 / T \ 

E(d w ) 


0.100 


0.100 


0.100 


0.100 


0.100 


0.100 


0.100 


0.100 


0.100 


0.100 


0.216 


0.172 


0.142 


0.118 


0.097 


0.080 


0.064 


0.050 


0.037 


0.024 


0.120 


0.120 


0.120 


0.120 


0.120 


0.080 


0.080 


0.080 


0.080 


0.080 


0.225 


0.176 


0.143 


0.117 


0.095 


0.077 


0.061 


0.047 


0.035 


0.023 


0.140 


0.140 


0.140 


0.140 


0.140 


0.060 


0.060 


0.060 


0.060 


0.060 


0.245 


0.189 


0.149 


0.116 


0.090 


0.069 


0.053 


0.040 


0.029 


0.019 


0.160 


0.160 


0.160 


0.160 


0.160 


0.040 


0.040 


0.040 


0.040 


0.040 


0.272 


0.206 


0.159 


0.121 


0.087 


0.053 


0.039 


0.029 


0.021 


0.014 


0.180 


0.180 


0.180 


0.180 


0.180 


0.020 


0.020 


0.020 


0.020 


0.020 


0.300 


0.226 


0.173 


0.129 


0.090 


0.028 


0.020 


0.015 


0.011 


0.007 


0.198 


0.198 


0.198 


0.198 


0.198 


0.002 


0.002 


0.002 


0.002 


0.002 


0.328 


0.245 


0.187 


0.139 


0.094 


0.003 


0.002 


0.002 


0.001 


0.001 


0.200 


0.089 


0.089 


0.089 


0.089 


0.089 


0.089 


0.089 


0.089 


0.089 


0.250 


0.171 


0.137 


0.112 


0.092 


0.075 


0.060 


0.046 


0.034 


0.023 


0.400 


0.067 


0.067 


0.067 


0.067 


0.067 


0.067 


0.067 


0.067 


0.067 


0.419 


0.134 


0.106 


0.087 


0.071 


0.058 


0.046 


0.036 


0.026 


0.018 


0.600 


0.044 


0.044 


0.044 


0.044 


0.044 


0.044 


0.044 


0.044 


0.044 


0.605 


0.091 


0.072 


0.059 


0.048 


0.039 


0.031 


0.024 


0.018 


0.012 


0.800 


0.022 


0.022 


0.022 


0.022 


0.022 


0.022 


0.022 


0.022 


0.022 


0.799 


0.046 


0.037 


0.030 


0.024 


0.020 


0.016 


0.012 


0.009 


0.006 


0.990 


0.001 


0.001 


0.001 


0.001 


0.001 


0.001 


0.001 


0.001 


0.001 


0.990 


0.002 


0.002 


0.002 


0.001 


0.001 


0.001 


0.001 


0.000 


0.000 
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p-dimensional sample vectors with a covariance matrix S. Suppose N > p. 
The sample covariance matrix 



N 



1 N 

i=i 



ar 



where a; is a sample mean vector, is an unbiased estimator of 17. We consider 
the estimation problem of r (defined by (pQ)) based on A. 
We define the following notations; 



v (* (jv) y y 



V (y (Ar) )' ) 



y 



a? 



(0 



a;, i = 1,...,JV, 



C^I^-iV -1 11', 



where Jjy is the iV-dimensional identity matrix, and 1 is the iV-dimensional 
vector with unit as each element. We find that 



CX 



The expression of A 



A 



1 



-Y'Y 



N-1 

is inconvenient, since the rows of Y are linearly constrained. Notice that 



A 



1 



N 



Y'Y 



1 



N-1 



-X'C 2 X 



N 



X'CX 



Using the decomposition of C = 0\0' x , 0\ G Vn-i,n, where Vn-i,n is the 
Stiefel Manifold, if we put Z = 0[X, then we have 



N-1 



Z'Z. 



(6) 



The distribution of A is determined by Z through (Q, where Z is not de- 
generated. 

The most frequently postulated situation is that Xi i — 1, . . . , N is inde- 
pendently distributed as the p-variate normal distribution N p (fi, U). Then 
(N — I) A is distributed as Wishart distribution W P (N — 1, U). This is dis- 
tiributionally equivalent to postulating that Z is distributed as 



n = N-l. 



The density function of Z is proportional to 

exp{-{l/2)tiZ'ZU- 1 )\U\~ n/2 . 



(7) 
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One of the natural generalizations of (J7J) is an elliptically contoured distri- 
bution, the density of which is given by 

fiiiZ'ZS- 1 )^- 11 ' 2 . (8) 

We formulate our estimation problem as follows; 
Z is a n x p (n > p > 2) random matrix, and its density with respect to 
Lebesgue measure on R np is given by (JBJ) with some function /(•) on R, where 
£ is an unknown positive definite p-dimensional matrix. We just observe 

S = Z'Z (9) 

We consider the estimation of the population contribution rates r given by 
(ED) based on S. 



2.2 Distribution of the Sample Contribution Rates 

From (jSJ) and (jHJ), the density of S is given by 

Cl /(tr5i:- 1 )|5| (n - p - 1)/2 |i:r n/2 (10) 

with some constant C\ (the proof can be found in Appendix). When fix) = 
exp(— x/2) this density function is that of a Wishart distribution. The dis- 
tribution S and the parameter £ are equivariant with respect to the trans- 
formations S — > BSB', £ — > B£B' for any p-dimensional nonsingular 
matrix B. 

The eigenvalues of S are denoted by k, 1 < i < p. We derive the 
distribution of the sample contribution rates 

A _ u 

ttj _p , I 1, . . . ,p. 

They are on the hyperplane 



V 



v 

|(d 1; ...,d p ) > • • • > d p > 0, = 1 1 . 

i=l 



We use the coordinate system d = (di, . . . , c? p -i) for P. The range of d is 
given by 

p-i 



XT 



|d di > ■ ■ ■ > d p -i > 0, y^dj < 1 j. 

i=l 



Theorem 1 Let 

S = HLH', L = diag(/ 1 ,...,/ P ), HeO{p) (11) 

be the spectral decomposition of S, where 0(p) is the set of p- dimensional 
orthogonal matrices. 
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(i) The density function of(d, H) with respect to the product measure between 
Lebesgue measure on RP' 1 and the invariant probability measure on 0(p) 
is given by 

c 2 \S\~ n/2 F{d) (trS^HDH'J , (12) 

where c 2 is a constant, D = diag(di, . . . , d p _i, 1 — £)?=i dj), and F(-) is a 
positive-valued function on T>* that is independent of /(•) in (jSJ) or S. 

(ii) The density function of d with respect to Lebesgue measure on RP~ l is 
given by 

c 2 f[r- n/2 F(d) [ (tiT- l HDH'\ nPl \ p ^dH), (13) 

~ = 1 JOM V J 



where T = diag(rx, . . . , r p ). 



Proof From (flQjl . the density function of I = (h, . . . , l p ) and H with respect 
to the product measure between Lebesgue measure on R p and on 0(p) 
is given by 

c Cl \S\- n ' 2 f[t- p ' l)/2 \[(k-h)f{^S- l HLH^ 

i=l i<j 

with some positive constant c. (See e.g. (22) on the pl05 of Muirhead(1982).) 
Let t = YTi=i h an d consider the transformation of coordinates 

I (t,d). 

Since the Jacorbian is given by J (I — > (t, d)) = t p ~ x y the density function of 
t, d and H with respect to dt x dd x [i p , p (dH) is given by 

c Cl \S\- n ' 2 t np/2 - l F(d)f(t txS' l H DH' ) , (14) 



where 



p-i p-i 



i<j<j<p— i i=i j=i 



X 

1=1 i=l 



Integrate (fl4"|) over {t | < t < oo}, then 

cd\S\- n/2 F(d) J t nr}/2 - l f(ta)dt (a = trS' 1 H DH' >0 

= cci|17|- n/2 F(d) a-" p/2 / t^^fitjdt. 

Jo 
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Now we have the density function of (d, H) with respect to dd x fi pp (dH) 

as 

c ci\E\- n/2 (j^ t np/2 -V(t)rf^F(d) (tiZ- 1 HDH'^ 



Substituting 



C Ci 



oo 

t np / 2 - l f{t)dt 







with C2, we have ( TTzj) . 

Integrating (fl2l) over 0(p), we have the density of d as 



r. 2 TTA- n/2 F(d) / (tiU- 1 HDH') nP/2 fi p>p (dH). (15) 
7=7 ■/o(p) v J 



Let 



2; = £r'Ai7, A = diag(A 1; ...,A p ), H e 0(p) 

be the spectral decomposition of S. Since /x PiP is the invariant probability on 
0(p), HH has the same distribution as H. Therefore (JTBT) equals 

c 2 TT A7 n/2 F(d) / ftrA- 1 ja r r>jy r, ) _np/ V PlP (rf-H'). 

7=1 ^o(p) v J 



Substituting A, with r, X]j=i ^j; we have f|T3|) . I 

It is noteworthy that the distribution of d is independent of /(•) in (JH1) 
and depends on £ only through the population contribution rates r. Note 
that Johnson and Grayvill (1972) deals with the distribution of the sample 
contribution rates when U = I p . 

2.3 New Estimator 

In order to derive a new estimator f(d) = (fi(d), . . . , r p (d)) of r = (r 1; . . . , r p ) 
that has a certain superiority to the classical estimator, we take a decision 
theoretic approach here, that is, we compare estimators via their risks with 
respect to a certain loss function. Straightforward approach is to use a loss 
function that directly measures the distance between r and its estimator f. 
However the exact property of the sample contribution rates d are difficult 
to derive under the assumption of small samples. Instead we evaluate the 
performance of r(d) as the components of an estimator of S. 

Consider an estimator of U combining r and the sample eigenvectors H 
in (ITT]) as follows; 

± = HT(d)H', f(d) = diag(fi(d), . . . , f p (d)). (16) 

Our approach is based on the following observations. According to Theorem 
[TJ the distribution of d is determined by r. Therefore we can suppose that 
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tiE = Y^h=i = 1 without loss of generality. In this case, the population 
contribution rates are equal to the population eigenvalues. In addition, the 
sample eigenvectors H are M.L.E. , hence consistent under the large sample 
asymptotics if /(•) in (10) is monotonically decreasing. (See Paul (2007) and 
Nadler (2008) for the discrepancy between the sample eigenvectors and the 
population counterparts for a large-dimensional matrix.) Therefore (Tl6|) is 
supposed to be a good estimator of S if r is a good estimator of r. 

The most common loss function about U and U is the entropy loss func- 
tion (Stein's loss function) 

L{±, S) = ti{SS- 1 ) - \og{\SS- l \) - p. (17) 

We evaluate the performance of £ through its risk with respect to this loss 
function. 

We consider one class of simple estimators given by 

f* = (t*,...,t;), t? = /3X i<i< P , (is) 

where (3*, i — 1, . . . ,p are positive constants. The classical estimator denoted 
by r° is given by 

r° = (ff,...,f;) s T? = 0?di, l<i<p, (19) 

where = 1, i — 1, . . . ,p. Correspondingly we define the two estimators 
S* and 27° as follows; 

S* = HT*H', T* = diag(f*, . . . , f*), (20) 
Z° = HT°H', f° = diag(f 1 °,...,7j). (21) 

We have the following result on the superiority of S* to S° . 

Theorem 2 If j3* (i = 1, . . . ,p) satisfy the following three conditions, then 
S* dominates S° with respect to the loss function ( [17]) . 
For some m (1 < m < p — 1), the next two inequalities hold; 

o <# <•••</?;< i < 0* m+1 < ■■■ < /?;, (22) 

m p 

J2{n + p-l-2i)(P*-l)+ (n + p + l-2i){/3*-l)<0. (23) 

i=l i=m-i-l 

Moreover, the third inequality 

E^T 1 < P (24) 

holds. 
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Proof Since both S* and S° is the function of d and H, their distributions 
are independent of /(•) from the result (i) of Theorem [TJ Therefore we can 
suppose f(x) = exp(— x/2), that is, S is distributed as a Wishart matrix; 

S~W p (n,E). (25) 

If S is distributed as in (1251) . the following Stein- Haff identity holds. (Ex- 
actly speaking, it is the application of Stein-Haff identity to an orthogonally 
equivariant estimator, see e.g. Lemma 2.1 of Dey and Srinivasan (1986)); 
Suppose S is decomposed as in ( ITT]) and S is given by 

S = Hdiag(Ml), <f> P (l))H', I = (h, ...,l p ). 

Then 

E[tr(tS- 1 )]=E[G{t,l)], 

where 

G(±l)=2 E M + 2 t^ + (»- P -l)tf. (26) 

l<j<j<p J 2=1 «=1 

If we use this identity, we have the following equation. 
E[L(E*,E)]-E[L(E°,E)] 

= E[ti{£*£- 1 ) - log E[ti{£° - log \S°\] (27) 
= E[G{±\ I) - G(S°, I) - log + log 
Substituting (ITS]) and (12"U]) into (12"61) . followed by simple calculation, we have 

^l^^J-v^ — r / _ / b — rZ^* - — / — 

^i=i J ' i<i<j< P H j ^i=i J ' i=i j 

1 p 

^"=1 •? i=l 

+ 2^/3;(l-4) + (n-p-l)^)9;} 
1=1 1=1 

^-0=1 •? i<«<j<p * J i=l 

+ 2E/3*(l-d i ) + (n-p-l)E/3*} 

i=l i=l 

P 



(28) 
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If we substitute 0* in ([28]) with 0? = 1 (1 < i < p), we have G(i; , Z). From 
these results, the inside of the brackets of the right-hand side in ( 1271) turns 
out to be 

(E^H 2 E (^-/3;) I -^ r + E(n + p + l-2z-2d,)(/3;-l)} 

i=l 

(29) 

Since /3* < /3? for 1 < i < j < p from (j22]) . (1291) is less than or equal to 

(e z 1 x> + p + 1 - 21 - 2d ^ (# - x ) - io s n ^- ( 3 °) 

(T22D says that 0* < 1 if 1 < i < m and that /3f > 1 if m + 1 < i < p. Using 
this fact together with the inequality < di < 1, Vz, we notice that (!30|) is 
less than or equal to 

P _^ m p p 

(E'i) {E( n +P- 1 - 2 ^)(^- 1 )+ E (n+P+l-2i)(0*-l)}-logl[0*. 

j=l i=l i=m+l i=l 

From (1231) . this is less than or equal to ^f=i Because of the 

inequality log(x + 1) < x, Vx > — 1, we have 

Emat 1 = Eiogi^;)- 1 - 1 + 1} < Ewr 1 - 1}, 

i=l i=l i=l 

which is nonpositive by (I24j) . I 



( 1221) of Theorem [2] means t* modifies the bias of the classical estimator 
which we mentioned in Section (TJ since lighter weight is given to di for the 
smaller z's and heavier weight for the larger z's. 

Choose an integer q such that 1 < q < p/2 — 1. Let 0^ (i — 1, . . . , p) be 
defined as 

n{n + p — 2q + 1 — 22)" 1 if 1 < 2 < m — q, 
0\ q) = { 1 if m-q+1 <i <p-m + q, (31) 

n(n + p + 2q + 1 — 2i) _1 if p — m + g + 1 < i < p, 

where m = [p/2], i.e., the largest integer that does not exceed p/2. Then 
0\ q) (i = l,...,p) satisfy (E2), ([23]) and (El}. In fact, (E2J) is clearly satisfied 
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from the definition. f l23|) and fl24|) are also satisfied as follows; 



Y j (n + p-l-2i)(B ( ? ) -1)+ £ (n+p+l-2^)(/3 l (9) -l) 



i=l 



i=m+l 



m—q 



= ^(n + p-l-2z)(/3 4 (9) -l)+ ^ (n + p + l-2z)($ ?) -l) 

j=l i=p—m+q+l 
m—q p 

<^(n + p-2q + l-2i)(B ( i q) + J] (n + p + 2g + 1 - 2z)(/3^ } - 1) 

i=l i=p— m+q+l 

m—q p 

= ^(-p + 2g-l + 2i) + Yl (-P-2g-l + 2i) 

i=l i=p— m+q+1 

= (m - + 2g - 1) + (m - g + 1) + (-p - 2g - 1) + (2p - m + q + 1)} 

= 0, 



i=l 

= Eio^r 1 - 1} 

L 

m—q p 

1 {^(p-2g + l-20+ £ ( p + 2g + l-2z)} 



8=1 



?7 



0. 



i=l 



j=p— m+q+l 



We give two examples of the estimators that satisfy the three conditions in 
Theorem [21 Let q=l, then 



(i) 



n(n + p — 1 — 2i)~ if 1 < i < m — 1, 

1 if m < i < p — m + 1, 

n(n + p + 3 — 2i) _1 if p — m + 2 < i < p, 



(32) 



The specific value of (3^ (1 < i < p) is given as follows; 
if p is even 









B {1) 


B (1) 


B {1) 

> J m+\ 


B (1) 

t J m+2 








n 


n 




n 
n+1 


1 


1 


n 
n— 1 




n 


n 


n+p— 3 


n+p— 5 


n— p+5 


n— p+3 



if p is odd 









Pm-l 


j->m 


Mm+l 


o(l) 


o(D 
Hm+3 








n 


n 




n 
n+2 


1 


1 


1 


n 
n-2 




n 


n 


n+p— 3 


n+p—5 


n— p+5 


n— p+3 
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The estimator 

£U = HTWH', f( 1 )=diag(f 1 (1) ,...,f«), (33) 

where = B^d^ 1 < i < p, dominates IJ° if p > 4. Another estimator 
that satisfies the three conditions of Theorem[2]is given by q = 2, which leads 
to 

{n(n + p-3-2i) _1 ifl<i<m-2, 
1 ifm-l<i<p-m + 2, (34) 

n{n + p + 5 — 2i)~ 1 if p — m + 3 < i < p. 

(2) 

The specific value of B\ (1 < i < p) is given as follows; 
if p is even 





pf 




Pm-2 


Pm-l 




B {2) 

> J m+2 


Pm+3 




4% 


pf 


n 


n 




n 
n+1 


1 




1 


n 
n-1 




n 


n 


n+p— 5 


n+p— 7 


n-p+7 


n—p+5 



if p is odd 





pf 






£ (2) 
Pm-l 




B {2) 

Pm+3 


/3 (2) 
^m+4 








n 


n 




n 
n+2 


1 




1 


n 
n-2 




n 


n 


n+p— 5 


n+p— 7 


n-p+7 


n—p+5 



Note that 

/?« < /?f ) < pt°\= 1) for 1 < i < m , 
pV> > Bf ] > Bf\= 1) for m + 1 < i < p. 

The estimator 

±m = H f(2) Hl f(2) = diag( +(2) ; . . . , ff), 

where tv = B^di, 1 < i < p, dominates 27° if p > 6. 

3 Simulation Study 

In this section, we examine by simulation other preferable properties of the 
new estimator f * = (f*, . . . , f*), 

t* = B-*d h l<i< P , 

where /3*'s satisfy §22$-§ZM, especially when B* = B^ or fif\ 1 < i < p. 
We use the notation = (ti \ . . . , Tp), 

T? ) = $ ) d i , 3 = 1,2, l<i<p. 



(35) 
(36) 
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Table 2: Risk w.r.t. Quadratic Loss 



n 


72 


T 3 


T 4 




76 


t- 7 


7- 8 


T9 


7-10 


~ 


~ 1 


* 2 


0.10 


o.io 


0.10 


0.10 


0.10 


o.io 


0.10 


0.10 


0.10 


0.10 


3.57 


2.11 


2.56 


0.11 


O.ll 


0.11 


0.11 


0.11 


0.09 


0.09 


0.09 


0.09 


0.09 


2.71 


1.57 


1.91 


0.12 


0.12 


0.12 


0.12 


0.12 


0.08 


0.08 


0.08 


0.08 


0.08 


2.20 


1.25 


1.53 


0.13 


0.13 


0.13 


0.13 


0.13 


0.07 


0.07 


0.07 


0.07 


0.07 


1.96 


1.13 


1.36 


0.14 


0.14 


0.14 


0.14 


0.14 


0.06 


0.06 


0.06 


0.06 


0.06 


1.85 


1.10 


1.31 


0.15 


0.15 


0.15 


0.15 


0.15 


0.05 


0.05 


0.05 


0.05 


0.05 


1.84 


1.14 


1.32 


0.16 


0.16 


0.16 


0.16 


0.16 


0.04 


0.04 


0.04 


0.04 


0.04 


1.87 


1.21 


1.38 


0.17 


0.17 


0.17 


0.17 


0.17 


0.03 


0.03 


0.03 


0.03 


0.03 


1.89 


1.27 


1.42 


0.18 


0.18 


0.18 


0.18 


0.18 


0.02 


0.02 


0.02 


0.02 


0.02 


1.91 


1.32 


1.46 


0.19 


0.19 


n i a 

0.19 


n i a 

0.19 


0.19 


0.01 


0.01 


0.01 


0.01 


0.01 


1.94 


1.37 


1.50 


0.20 


0.09 


0.09 


0.09 


0.09 


0.09 


0.09 


0.09 


0.09 


0.09 


2.71 


1.78 


2.15 


0.30 


0.08 


0.08 


0.08 


0.08 


0.08 


0.08 


0.08 


0.08 


0.08 


2.86 


1.92 


2.31 


0.40 


0.07 


0.07 


0.07 


0.07 


0.07 


0.07 


0.07 


0.07 


0.07 


2.97 


2.02 


2.42 


0.50 


0.06 


0.06 


0.06 


0.06 


0.06 


0.06 


0.06 


0.06 


0.06 


3.09 


2.12 


2.53 


0.60 


0.04 


0.04 


0.04 


0.04 


0.04 


0.04 


0.04 


0.04 


0.04 


3.24 


2.23 


2.66 


0.70 


0.03 


0.03 


0.03 


0.03 


0.03 


0.03 


0.03 


0.03 


0.03 


3.44 


2.39 


2.84 


0.80 


0.02 


0.02 


0.02 


0.02 


0.02 


0.02 


0.02 


0.02 


0.02 


3.65 


2.56 


3.03 


0.90 


0.01 


0.01 


0.01 


0.01 


0.01 


0.01 


0.01 


0.01 


0.01 


3.94 


2.79 


3.28 



3.1 Risk Comparison 

We will compare the estimator and with the classical estimator 
through their risks with respect to the quadratic loss function; 

QL(f,r) = ^(l-f 4 /r 4 ) 2 . (37) 
i=i 

According to Theorem [21 the plug-in estimator made from j = 

1.2 dominates another plug-in estimator IJ^ from with respect to the 
entropy loss function. We are interested in a more direct comparison among 
f(°), f « and f ( 2 ) using ((37]). 

We generated 10000 random 10-dimensional Wishart matrices with the 
degree of freedom 30 under several patterns of the population contribution 
rates, t — (r l5 . . . ,r 10 ). The Table [2] shows the simulation result, where the 
first 10 numbers in each row are the population contribution rates and the last 
three numbers are the simulated risks for the three estimators j = 0, 1, 2 
(all the numbers are rounded to the second decimal place). The risk of 
is smaller than that of f^ ^ by 30% to 40%. Since is located between 
-f(°) and (see f|35]) ). its risk reduction is smaller than Nevertheless 
it still reduces the risk by 17% to 30% compared to r^. From these results, 
we can conclude that the new estimators are substantially improved over the 
classical estimator in view of the quadratic risk. 
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3.2 Estimation of Dimension 



As we mentioned in Section CD, r, the population contribution rates, is one 
of the most basic tools for deciding the dimension in principal component 
analysis or factor analysis. As the first step in deciding the dimension, the 
choice of an estimator r for r is an important task, hence we are interested 
in how the new estimator, t^*' , makes a difference compared to the classical 
estimator, t^°\ in the decision of the dimension. 

Suppose that generated in the following m-factor 

model ; 

x = a + Bz + e, (38) 

where a is a constant p-dimensional vector, B is a p x m (p > m) factor 
loading constant matrix with the rank of m, z is a m-dimensional random 
factor, and e is the p-dimensional error term which is independent of z. If 
we suppose the covariance matrices of z and e are respectively given by 



V(z) = 27 , V(e) = a 2 l 



then the covariance matrix S of x equals 

27 = BS G B' + a 2 I p . (39) 

If we denote the eigenvalues of BS B' by (1 <i < p), then (1 < i < p), 
the eigenvalues of U, are given by 



& + cr 2 if i = 1, . . . ,m, 

2 . f ■ , , ( 4 °) 

a if i = m + 1, . . . ,p 



since £j = 0, (m+1 < i < p). In the large- dimensional random matrix theory, 
the model (I3"5|) is called "spiked covariance model" , since the part of & + a 2 
seems spiked into the long flat part cr 2 . If z and e have normal distributions, 
then x is also normally distributed with the covariance matrix (1391) . hence 
we can suppose that S ~ W p (n, £). 

We made a simulation under the condition 



S~W w (n,A), n = 30, 100, A = diag(A x , . . . , A p ), 

where Aj (1 < i < p) is given by (HOI with ten patterns of (m, £i, . . . ,£ m ) 
and a 2 fixed to be unit. S is generated 10000 times, and for each time 
we recorded the dimension decided from the six methods composed by the 
combination of two criterions in SectionT] and three estimators t^°\ 
t 1 - 2 ); the two criterions are "the cumulative percentage of the (estimated) 
population eigenvalues" with t* = 0.8 in (j2j) (say criterion 1) and "the relative 
size of each (estimated) population eigenvalue" (say criterion 2). 

Table [3] and H] are the result of the simulation. To explain the meaning 
of each number in the table, take the first case in Table [3] as an example, 
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Table 3: Histogram of Estimated Dimension, n=30 



Case 1 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


19 


19 


19 


19 


19 


1 


1 


1 


1 


1 


c.i k fw 








104 


9895 


1 

















C.l k fW 











9 


9991 

















C.I k 











339 


9661 

















C.2 & fW 








285 


6480 


3235 

















C.2 & f W 








420 


6346 


3234 

















C.2 & f ^ 








2275 


4930 


2795 

















Case 2 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


23.5 


23.5 


23.5 


23.5 


1 


1 


1 


1 


1 


1 


c.i & fw 








7348 


2652 




















C.l & f w 











9828 


172 

















C.l & 








17 


9982 


1 

















C.2 & fW 





1 


1087 


8912 




















C.2 & fW 





5 


1348 


8647 




















C.2 & f ^ 





2 


3801 


6197 




















Case 3 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


31 


31 


31 


1 


1 


1 


1 


1 


1 


1 


c.i k fw 





376 


9624 























C.l & 








2410 


6641 


948 


1 














C.l & 








9987 


13 




















C.2 k fW 





22 


9978 























C.2 & 





50 


9950 























C.2 k 





30 


9970 























Case 4 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


46 


46 


1 


1 


1 


1 


1 


1 


1 


1 


c.i k fw 





10000 


























C.l & 





1 


137 


3545 


5907 


410 














C.l & f ^ 





8008 


1903 


88 


1 

















C.2 & fW 





10000 


























C.2 k f W 





10000 


























C.2 & 





10000 


























Case 5 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


91 


1 


1 


1 


1 


1 


1 


1 


1 


1 


c.i k fw 


9991 


9 


























C.l & fW 











2 


776 


8170 


1052 











C.l & f ^ 


1184 


3743 


3916 


1064 


93 

















C.2 & 


10000 





























C.2 & 


10000 





























C.2 k 


10000 
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Case 6 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


27 


23 


19 


15 


11 


1 


1 


1 


1 


1 


C.l k f W 








1018 


8981 


1 

















c.i k f w 











48 


9952 

















C.l k 











1535 


8465 

















C.2 & 





2 


1483 


7619 


896 

















C.2 & fW 





12 


1922 


7171 


895 

















C.2 k f ^ 





2 


5407 


3904 


687 

















Case 7 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


32 


25.5 


19 


12.5 


6 


1 


1 


1 


1 


1 


C.l & fW 








5931 


4069 




















c.i k fw 











688 


9311 


1 














C.l & f ^ 








2 


6764 


3234 

















C.2 k fW 





35 


4702 


5254 


9 

















C.2 & fW 





107 


5209 


4676 


8 

















C.2 & f ^ 





49 


8246 


1701 


4 

















Case 8 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


37 


28 


19 


10 


1 


1 


1 


1 


1 


1 


C.l k t w 





36 


9847 


117 




















c.i & fw 











9013 


987 

















C.l k 








910 


9090 




















C.2 & 





107 


7496 


2397 




















C.2 & fW 





317 


7705 


1978 




















C.2 & f ^ 





168 


9347 


485 




















Case 9 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


46 


31 


16 


1 


1 


1 


1 


1 


1 


1 


C.l & f w 





4471 


5529 























c.i k f w 








598 


6457 


2923 


22 














C.l & f ^ 








9971 


29 




















C.2 & f W 





1022 


8978 























C.2 & fW 





1756 


8244 























C.2 & f ^ 





1234 


8766 























Case 10 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, • • • , Aio 


81 


11 


1 


1 


1 


1 


1 


1 


1 


1 


C.l k f w 


6232 


3767 


1 























c.i & fw 











14 


3144 


6766 


76 











C.l & f ^ 





4482 


4664 


824 


30 

















C.2 k fW 


3939 


6061 


























C.2 & 


5952 


4048 


























C.2 k 


5171 


4829 
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Table 4: Histogram of Estimated Case, n=100 



Case 1 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai ..... Ai n 

15 1 _1_U 


19 


19 


19 


19 


19 


1 


1 


1 


1 


1 


C.l & f w 











9536 


464 

















C.l & f W 











1259 


8741 

















C.l & f ^ 











1922 


8078 

















C.2 & fW 











50 


9950 

















C.2 & 











50 


9950 

















C.2 & 











50 


9950 

















Case 2 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai ..... Ai n 

15 7 _1_U 


23.5 


23.5 


23.5 


23.5 


1 


1 


1 


1 


1 


1 


C.l & 








287 


9713 




















C.l & f ^ 











10000 




















C.l & 








9 


9991 




















C.2 & f(°) 








1 


9999 




















C.2 & f ^ 








1 


9999 




















C.2 & f ^ 








2 


9998 




















Case 3 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 




31 


31 


31 


1 


1 


1 


1 


1 


1 


1 


C.l & 








10000 























C.l & 








10000 























C.l & f ^ 








10000 























C.2 & 








10000 























C.2 & 








10000 























C.2 & f ^ 








10000 























Case 4 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai ..... Am 


46 


46 


1 


1 


1 


1 


1 


1 


1 


1 


C.l & 





10000 


























C.l & 





10000 


























C.l & 





10000 


























C.2 & 





10000 


























C.2 & f ^ 





10000 


























C.2 & f ^ 





10000 


























Case 5 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


91 


1 


1 


1 


1 


1 


1 


1 


1 


1 


c.i & f w 


10000 





























C.l & f w 


9995 


5 


























C.l & f ^ 


9999 


1 


























C.2 & 


10000 





























C.2 & f W 


10000 





























C.2 & 


10000 
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Case 6 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


27 


23 


19 


15 


11 


1 


1 


1 


1 


1 


c.i & f w 











9998 


2 

















C.l & f W 











8849 


1151 

















C.I & f ^ 











9422 


578 

















C.2 & fW 








15 


5304 


4681 

















C.2 & fW 








20 


5299 


4681 

















C.2 & f W 








87 


5239 


4674 

















Case 7 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


32 


25.5 


19 


12.5 


6 


1 


1 


1 


1 


1 


C.l & fW 








1113 


8887 




















C.l & f W 








1 


9999 




















C.l & 








25 


9975 




















C.2 & fW 








1095 


8905 




















C.2 & f W 








1228 


8772 




















C.2 & 








2321 


7679 




















Case 8 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . • , Aio 


37 


28 


19 


10 


1 


1 


1 


1 


1 


1 


C.l & f w 








9964 


36 




















C.l & f w 








5754 


4246 




















C.l & 








8784 


1216 




















C.2 & 








6280 


3720 




















C.2 & f W 








6531 


3469 




















C.2 & f ^ 








7859 


2141 




















Case 9 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . • , Aio 


46 


31 


16 


1 


1 


1 


1 


1 


1 


1 


C.l & fW 





1628 


8372 























C.l & f w 








10000 























C.l & 





49 


9951 























C.2 & 





19 


9981 























C.2 & f W 





38 


9962 























C.2 & 





22 


9978 























Case 10 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Ai, . . . , Aio 


81 


11 


1 


1 


1 


1 


1 


1 


1 


1 


C.l & fW 


6679 


3321 


























C.l & fW 


208 


9792 


























C.l & 


1108 


8892 


























C.2 & fW 


2952 


7048 


























C.2 & fW 


4039 


5961 


























C.2 & 


3612 


6388 
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where A*, z = 1, . . . , 10 are given by 19, 19, 19, 19, 19, 1, 1, 1, 1, 1. (Note that 
for every case, Aj's are designed so that X^=i ^* = hence each A« equals 
Ti in percent figures.) The row of "C.I & r^" is the histogram of the 
estimated dimension by the combination of criterion 1 and . The boldface 
position ( in this case "1" ) indicates that the "true" dimension decided by 
criterion 1 from the population eigenvalues, Aj's. In this simulation, the 
"true" dimensions are designed to take the same value by either criterion 1 
or 2. 

We can observe following points from Table [3j With respect to the cri- 
terion 1, we notice that the classical estimator, f(°\ tends to underestimate 
the dimension (see Case 1, 2, 7, 9), while tend to overestimate it (see 
Case 3-10). From (135]) . we notice that is located between and f^K 
Though is still likely to overestimate the dimension (see Case 5, 6, 8, 10), 
the tendency is weakened compared to . On the criterion 2, we can not 
find as significant a difference as criterion 1 among three estimators. Every 
estimator tends to underestimate the dimension in some cases (see Case 1, 
7) and overestimate it in other cases (see Case 10). 

In most cases in Table HJ the estimation for dimension is made correctly. 
However, despite a high degree of freedom, we still observe the tendency of 
t(°) to underestimation (see Case 1) and that of or to overestimation 
(see Case 10) with respect to the criterion 1. 

Though both underestimation and overestimation are undesirable, the 
former is more crucial, and have more substantial effect on the results ob- 
tained in principal component analysis or factor analysis, since important 
component (factor) is neglected (see e.g. the comment in p 278 in Fabriger 
et. al. (1999)). In this sense, and are superior to . The tendency 
to overestimation of is weakened in since /?f^'s are closer to unit. 
We can correct the overestimation further by selecting that are much 
closer to unit, but still satisfy the three conditions in Theorem [2j 

4 Conclusion 

We can summarize the results of this paper as follows; 

1. The distribution of the sample contribution rates is identical within a 
family of elliptically contoured distributions. It is determined solely by 
the population contribution rates. 

2. A class of new estimators of the population contribution rates was 
derived. In the estimation of the normalized population covariance 
matrix, the estimator composed of the new estimator and the sample 
eigenvectors dominates the estimator composed of the classical esti- 
mator and the sample eigenvectors with respect to the entropy loss 
function. 
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3. A simulation study shows that the new estimators perform substan- 
tially better than the classical estimator with respect to the risk de- 
rived from a quadratic loss function. Another simulation study shows 
that the new estimators tend to overestimate the dimension. They are 
more suitable than the classical estimator for the decision of dimension 
in principal component analysis or factor analysis, since the classical 
estimator is likely to underestimate the dimension. 
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6 Appendix 
6.1 Proof of © 

We can suppose S ~ W p (n, A), A = diag(Ai, . . . , X p ), since the distribution 
of (ij's depends only on A. Let A = (5ij) be defined as 



(3) of Lawley (1959) gives the following expansion of /* (1 < r < p), 



Since Ef=i ^* = Ei=i(Aj + &a)i we have the Taylor expansion of (E?=i ^t) 1 



A = n S - A = A- A. 




) 



+ 0(||<5|| 4 ). 



(41) 



1 1 



ELii? Ef=i(Ai + da) 



Ej<j figfijj 



(42) 



EL: si E 



l<i,j,k<p °ii jj°kk 



+ o(\\s\\ 4 )- 



(ELi a,) 3 



(Eti a,) 4 
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Combining (j4ip and (142]) . we have 



r v^p 7* 
Z-/i=l l i 



A r 1 A r 



(43) 



We can easily calculate the low-dimensional moments of 5. They are given 
as follows; 

£N =0, 1 < z < p. (44) 

™= o, if^J: (45) 



= n^W, \H±j. (46) 



In addition, 



E[\\8\\ k ] = 0{n- {k - 1] ). (47) 
Substituting (l4"4"|) - ( 1471) into ( T431) . we have the desired result. I 

6.2 Proof of ([TDD 

Let the singular value decomposition of Z denoted by 

Z = ODH, D = diag(di, . . . , dp), di > ■ ■ ■ > d p > 0, O 6 V^, n , e 0(p), 

where O(p) is the set of p-dimensional orthogonal matrices. From the Jaco- 
bian of this decomposition (see e.g. Theorem 5 of Uhlig (1994)), we have 

v 

dZ oc Y[ d n ~ p J](d 2 - df) n P A dH ) V P ,n{dO) dd, 

i=l i<j 

where and /i Pi „ are the invariant probability measures respectively on 
0(p) and V p>n . Further by the transformation di — > U = d\ i = 1, . . . ,p, we 
have 

v 

dZ oc \[t { r p - 1)/2 \[{U - tj) li PtP (dH) fi P , n (dO) dt. 

i=l i<j 

Notice that 

S = H'TH, T = diag(ti ,...,t p ), (48) 
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hence, 

/(trZ'Z2r 1 )|27|- n / 2 dZ 

oc /(trSIT 1 )|S| (n - p - 1)/2 |i;|- ,l/2 Y[{ti -tj) n PlP (dH) fi Pt „{dO) dt. 

i<j 

If we integrate the right side of the above equation over V p>n , we have the 
following density function of H and t with respect to n P>p {dH)dt 

c /(trSir 1 )|S|( n - p - 1 )/ 2 |i;|- n / 2 Y[(ti - tj) (49) 

with some constant cq. Combined with a formula about the spectral decom- 
position (j48D (see e.g. (22) on the pl05 of Muirhead(1982)) 

dS oc Y\_(Ji — tj) ti> PP (dH) dt, 

i<j 

fj49|) leads to the following density of S with respect to Lebesgue measure 

c l f{tTSZ- 1 )\S\ {n -' p - 1)/2 \Z\~ n/2 (50) 
with a normalizing constant C\. I 
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